
function [Z, res] = ALSbaseline(Y, lambda, p, itter)
% Estimate baseline with asymmetric least squares
% 
% Reference:
% Eilers, P. and H. Boelens "Baseline Correction with Asymmetric Least
% Squares Smoothing."
% http://www.science.uva.nl/~hboelens/publications/draftpub/
%
% INPUT:    Y       = Signal
%           lambda  = smoothness     
%                           "and 10^2 > lambda > 10^9, but exceptions
%                               may occur." 
%           p       = asymmetry      
%                           "We found that generally 0.001 > p > 0.1 is a
%                               good choice (for a signal with positive
%                               peaks)
%           itter   = number of interratiosn    (DEFAULT = 10)
%
% OUTPUT:   z       = calculated baseline
%           res     = residuals (deviation of baseline from signal)

% (c) F. Dworkowski
% v2.0, 2012-02-21 : now supports multiple spectra in one input

if nargin<4,
    itter = 20;
    if nargin<3
        p = 0.01;
        if nargin<2
            lambda = 10e5;
            if nargin<1
                help ALSbaseline
                return
            end
        end
    end
end

if itter < 1 || itter > 1000
    ii = 10;
else ii = itter;
end
    
[m,n] = size(Y);
D = diff(speye(m), 2);
for i = 1:n
    w = ones(m, 1);
    y=Y(:,i);
    for it = 1:ii
        W = spdiags(w, 0, m, m);
        C = chol(W + lambda * D' * D);
        z = C \ (C' \ (w .* y));
        w = p * (y > z) + (1 - p) * (y < z);
    end
    Z(:,i)=z; 
end
res = Y - Z;
